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A mathematical model was developed to simulate flow control applications using plasma 
actuators. The effects of the plasma actuators on the external flow are incorporated into 
Navier Stokes computations as a body force vector. In order to compute this body force 
vector, the model solves two additional equations: one for the electric field due to the applied 
AC voltage at the electrodes and the other for the charge density representing the ionized 
air. The model is calibrated against an experiment having plasma-driven flow in a quiescent 
environment and is then applied to simulate a low- pressure turbine flow with large flow 
separation. The effects of the plasma actuator on control of flow separation are 
demonstrated numerically. 
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Nomenclature 

= pressure coefficient 
= axial chord for Pak-B blade 
= elementary charge, C 
= electric field, N/C 
= permittivity, e = s r s o 
= relative permittivity 

= permittivity of free space, 8.854xl0' 12 C I 2 /Nm 2 
= body force vector, N/m 3 
= freestream turbulence intensity, % 

= total electric potential, Volt, (& = <j) + (p 
= electric potential due to external electric field. Volt 
= electric potential due to net charge density, Volt 
= Boltzmann’s constant 
= length of the electrode 
= Debye length, m 

= location parameter for Gaussian distribution 
= ion density in the plasma 
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n e 

— 

electron density in the plasma 

n o 

= 

background plasma density 

CO 

= 

frequency, Hz 

Re 

= 

Reynolds number based on inlet velocity and axial chord 

Pc 

= 

net charge density, C/m 3 

a 

= 

scale parameter for Gaussian distribution 

t 

= 

time, s 

T 

= 

temperature of the species 

u in 

= 

inlet velocity 

u 

= 

streamwise velocity 

x.y 

= 

coordinates 

Yn 

= 

wall normal distance 


D ec 
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I. Introduction 

ecent advances in the area of flow control have led to the application of new flow control strategies to various 


ow fields. 3 One of the major areas for applications of flow control is the gas turbine, which has the potential 
to lead to significant benefits in improved design and performance. 4 ' 3 Several researchers have investigated the 
effectiveness and feasibility of different flow control methods for reducing losses and improving performance in gas 
turbine engine component flowfields. 6 " 9 Although many of these flow control methods have been beneficial in 
laboratory tests, the incorporation of such systems into actual applications is not straightforward. It involves 
development and testing of the control system for the actuators, as well as the hardware development of the device. 
An experimental approach requires numerous costly and time consuming trial-and-error iterations. With advances in 
computational fluid dynamics (CFD) and computer hardware, CFD offers an alternative method for designing and 
understanding complicated flow control systems. Continued development of CFD simulation tools and control 
algorithms is essential in order to adapt the available flow control technologies to design of modem gas turbine 
engines. 

Among the active flow control techniques the use of plasma actuators has been demonstrated to be effective in 
several applications including flow separation and boundary layer control. The plasma actuators consist of two 
electrodes that are located on surface separated by a dielectric material as shown in Fig. 1. A high-voltage AC 
supplied to the electrodes causes the air in their 
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Figure 1. Schematic of plasma actuator. 


vicinity to weakly ionize. The ionized air (plasma) 
in the presence of the electric field gradient 
produced by the electrodes, result in a body force 
vector acting on the external flow that can induce 
steady or unsteady velocity components. The 
effectiveness of plasma actuators in controlling flow 
separation has been demonstrated by several 
researchers. Examples of various applications 
include boundary layer control by Roth et al. 10 , 
exciting boundary layer instabilities on a sharp 
cone at Mach 3.5 by Corke et al. 11 , lift augmentation 
on a wing section by Corke et al. 12 , wing section leading edge separation control by Post and Corke 13 , unsteady 
vortex generation by Orlov et al. 14 , separation control on low pressure turbine blades by Hultgren and Ashpis 13 , 
Huang et al. 9 and Jacob et al. 16 , separation control on stationary and oscillating airfoils by Post and Corke 17 , and tip 
clearance control in turbines. 18 The results from these experiments indicate that several parameters have to be taken 
into consideration for effective flow control. The important parameters are location of the actuators on the surface, 
orientation, size, and relative placement of the embedded and exposed electrodes, applied voltage, and frequency of 
the actuation. The many parameters affecting the performance of the actuators makes incorporation of such systems 
into actual applications and the development and testing of optimized control systems for actuators a complicated 
task. CFD simulation tools can enable design and optimization of such complex flow control systems. 
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Computational studies of plasma flow control have been limited in comparison to the vast number of 
experimental studies. Examples of computational work include modeling of the effects of the plasma actuator on the 
external flow as a time averaged mean body force distribution in a triangular area located above the embedded 
electrode by Shyy et al. !9 , modeling the effects of plasma actuators using a potential flow method based on a doublet 
model by Hall et al. 20 , and simulation of plasma actuator effects using a lumped parameter model by Orlov and 
Corke 21 . In contrast to these simpler approaches, Roy and Gaitonde 22 23 modeled the collisional surface dielectric 
barrier discharge at atmospheric pressures by solving the charge continuity, charge momentum, electric field and 
potential distribution. In their approach helium was used as the carrier gas because its plasma coefficients and 
chemistry were available and the electric force obtained as a product of the charge and electric field was used to 
calculate the flow actuation effect. 

In this paper, we outline a new numerical simulation methodology for active flow separation control applications 
using plasma actuators. This is computationally efficient and easy to implement into existing RANS codes. In this 
new approach, the effect of the plasma actuators on the external flow is incorporated into Navier Stokes 
computations as a body force vector. The body force is obtained as a product of the net charge density and the 
electric field. The new model solves the Maxwell equation to obtain the electric field due to the applied AC voltage 
at the electrodes and an additional equation for the charge density representing the plasma density. Details of the 
model development and implementation are provided below. Results of computations for a single plasma actuator 
in quiescent flow and computations of flow separation control over a Pak-B low-pressure turbine blade using a 
single plasma actuator are given in section 3. Concluding remarks are provided in section 4. 


II. Plasma Physics Modeling and CFD Computations 

The plasma actuators are formed by a pair of electrodes separated by a dielectric material. The actuator is placed 
in the surface with one electrode exposed to the surroundings and the other one embedded in the surface below the 
dielectric material (Figure 1). When a high-voltage AC is supplied to the electrodes, this arrangement causes the air 
in their vicinity to weakly ionize. The ionized air in the presence of the electric field gradient produced by the 
electrodes, results in a body force vector acting on the external flow that can induce steady or unsteady velocity 
components. This body force can be expressed in terms of the applied voltage and incorporated into the Navier 
Stokes equations. By neglecting magnetic forces, the electrohydrodynamic (EHD) force can be expressed as 

fs=P c E (D 


where, f g is the body force per unit volume, p c is net the charge density and E is the electric field. If the time 
variation of the magnetic field is negligible, as is often the case in plasma, the Maxwell’s equations give rise to 
Vx E ~ 0. This implies that the electric field can be derived from the gradient of a scalar potential 10 ' 24 : 


E = (2) 

Gauss’s law yields: 

V ■{£&) = p c (3) 


or 

V-(eV<t>) = -p c (4) 

where e is the permittivity. The permittivity can be expressed as: 

S = £,.£„ (5) 

where e r is the relative permittivity of the medium, and e o is the permittivity of free space. 


There still remains some controversy as to the physical mechanism by which the plasma actuator modifies the 
surrounding flow field. While it is clear that the body force developed via equation (1) has a large impact on the 
flow field, there is still some debate as to the effect of local temperature rise and reduced viscosity on the boundary 
layer. This accounts for the increase in velocity in the plasma region by increasing the shear rate to keep the viscous 
stress constant, as suggested by Leonov et al. 25 . Recent measurements appear to discount this reasoning, due to the 
large impact the small plasma region has on the surrounding flow 16 . Instead, the plasma region acts much like a 
source term in this regard. 
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Since the gas particles are weakly ionized, we can assume the potential <t> can be decoupled into two parts: one 
being a potential due to the external electric field, <p , and the other being a potential due to the net charge density in 

the plasma, (p , 

<& = Q + (p. (6) 

This approach is similar to the one used in numerical simulation of electroosmotic flows in which case the external 
electric field generates a force on the charged particles creating flow 26 . 

If we assume that the Debye thickness is small and the charge on the wall is not large, the distribution of charged 
species in the domain is governed by the potential on the wall due to the electric charge on the wall and is affected 
very little by the external electric field. Therefore, we can write two separate equations in terms of these two 
potentials, one for the external electric field due to 
the applied voltage at the electrodes: 

V-(£ r V0) = O (7) 

and another one for the potential due to the 
charged particles: 

V-(£,.Vp) = -(p c /£„). (8) 

Equation (7) is solved for the electric 
potential, 0, using the applied voltage on the 
electrodes as boundary conditions. The boundary 
conditions and the computational domain for 
Equation (7) are shown in Figure 2 for a single 
pair of electrodes. Equation (7) is solved with the 
appropriate e r value on both the air side and the 
wall. For air, £,. =1.0. We chose Kapton as the 
dielectric material, as in the experiments and it 
has an £ r value of 2.7. On the wall-air interface 
harmonic mean of £,> and£ r , must be used in order to conserve electric field as shown in Figure 3. In this figure 
point 1 is on air side and point 2 is in the dielectric material, with point / on the interface. 


Air Side: On Outer Boundaries: 

8<P . 

v(* rl v*)=o aT° 

s rl = 1.0 



Figure 2. Computational domain and boundary 
conditions for Equation 7. 
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Figure 3. The interfacial permittivity. 


At the dielectric wall and air interface, the electric field is conserved as shown in figure 3 by writing: 


f-1 =-.l 
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( 9 ) 
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or: 


d>.-& <P,-<Pi 



( 10 ) 


By using the electric analogy as shown in figure 3, the interfacial permittivity, e f , should take the harmonic mean 
of the permittivity of air, £, , and the permittivity of dielectric material, £, : 
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Applied AC voltage imposed at the exposed (upper) electrode as boundary condition is: 

m=r x f{t) 


( 12 ) 


The wave form function, / ( t ), can be a sine wave given by: 

/(/) = sin(2/ra>0 (13a) 

or it can be a square wave given by: 

,1 for sm{2kox)>Q 

/W={_i for sin(27Tfitf)<0 ( jb) 

where co is the frequency and is the amplitude. The embedded electrode is prescribed as ground by setting the 
electric potential to zero on that electrode. At the outer boundaries, 3<p/dn = 0 is assumed. 


Next we consider the electric potential due to the net charge in the plasma, equation (8). The net charge density 
within the plasma at any point is 24 : 


p c / £ 0 =e(n, -n e )l£ 0 ~ -^ 2 -^snp[e(p! kT t ) + exp(e<p/ kT e )] 


(14) 


where, e is elementary charge; and n e are the ion and electron densities in the plasma; n a is the background 
plasma density; k is the Boltzmann’s constant; T is the temperature of the species; and <p is the local electric 
potential. Expanding the exponential functions in a Taylor series for^)« T , equation (14) becomes, to lowest order 
of <p/T \ 


p c ! e o =-(e 2 nJ 0[(l/ «1) + (1/ kT')]<p. (15) 

Introducing the Debye length, , 


{\l^) = {e\l£ 0 )[[\lkT i ) + {\lkT ! )\ 

(16) 

equation (15) becomes; 


pje 0 =(-l/A l)q> 

(H) 

And, alternatively, the electric potential is: 


p = (-pA 2/ £„)- 

(18) 


Substituting (18) in Equation (8) we obtain an equation in terms of the net charge density at any point within the 
plasma. 

V(£V A .) = p ( ./Aj. (19) 
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Equation (19) is solved for the net charge 
density, p , , only on the air side of the domain 
using the boundary conditions shown in Figure 
4. A zero normal gradient for the net charge 
density is imposed on the solid walls except the 
region covering the lower electrode and the 
charge density is set to zero on the outer 
boundaries. Downstream of the exposed 
electrode, on the region over the embedded 
electrode the charge density is prescribed in 
such a way that it is synchronized with the time 
variation of the applied voltage <p{t) on the 
exposed electrode, /(/) in equation (12): 
p c Jx,t) = p™G(x)f(t) (20) 

where />” ax is the maximum value of the charge 
density allowed in the domain (in Coulomb/m 3 ), 
a parameter to be determined later. The variation of the charge density on the wall, p c lv , in the streamwise direction, 

x, is prescribed by a function G(x) which can be chosen to resemble the plasma distribution over the embedded 
electrode. Experimental results 24,27 28 suggest that this distribution is similar to a half Gaussian distribution given by 

G(x) = exp^-(,r-/r) 2 /(2er 2 )J (21) 

for x > 0 . In equation (21) p is the location parameter indicating the x location of the maximum, and a is a scale 
parameter determining the rate of decay. Choosing a very large a value results in a constant p c w distribution. The 
variation of G(x) with various values of a defined in terms of the length of the electrode, L e , is illustrated in figure 
5. In this figure p is equal to zero, and as a increases the decay rate decreases. 

In order to solve Eqs. (7) and (19), the Debye length, 

A, , maximum charge density on the wall, p ™** , scale 
parameter for charge density distribution, cr, have to be 
prescribed. These parameters can be calibrated based on 
a simple experiment having a pure plasma-driven flow in 
a quiescent environment, as will be explained in the next 
section. In the computations we have chosen the location 
parameter// such that the peak corresponds to the left 
edge of the embedded electrode as shown in Figure 4. 

Moreover, we have assumed that a takes a value of 0.3 
to allow a gradual decay of the charge density 
distribution from the left edge to the right edge. The 
values of the frequency and amplitude of the applied 
voltage, <D and <t> mm in equations (12), (I3a) and (13b) 
are given from the experiments. 

Once <p and p c are obtained from solutions of 
equations (7) and (19) respectively, the resulting body 
force vector is computed from Equation ( 1 ), that is: 
f B =p c E = p c (-V<p). (22) 

It must be noted that solution of equations (7) and (19) is done only once at the beginning of computations, since 
these equations do not contain a time derivative term. This can be achieved by normalizing <p using the value of AC 



Figure 5. Half Gaussian distribution given by 
Equation (21) 



Figure 4. Computational domain and boundary conditions 
for Equation 19. 
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voltage of the exposed electrode, <f >(/) , in equation (12). Equation (7) can then be solved by imposing a constant 
boundary condition equal to unity at the upper electrode. Once the dimensionless <j> distribution is determined, the 
dimensional 0 values at any given time can be obtained by multiplying this distribution with the corresponding 
value of <j)(t) given by equation (12). Similarly, equation (19) can be solved at the beginning of computations by 
solving the dimensionless charge density distribution, which is normalized by p™ a f{t) ■ This implies that the 
boundary condition for the dimensionless charge density on the wall region covering the embedded electrode is 
G(x). Once the solution for the dimensionless charge density is established, the dimensional values at any time can 
be calculated from this distribution by multiplying it with the corresponding value of p™*f(t) ■ 

This modeling approach for plasma actuators is implemented in the GHOST code developed at University of 
Kentucky. GHOST is a pressure-based code based on SIMPLE algorithm with second order accuracy in both time 
and space. Advection terms are approximated by a QUICK scheme and central differencing is used for the viscous 
terms. The ‘Rhie and Chow’ momentum interpolation method 29 is employed to avoid checkerboard oscillations 
usually associated with the non-staggered grid arrangement. This code is capable of handling complex geometries, 
moving and overset grids, and includes multiprocessor computation capability using MPI. Since multiple processors 
are used during the computations, it is more efficient to divide the computational domain into several smaller pieces 
with very fine grids and distribute the zones to processors with the consideration of equaling load balancing. The 
overset grid capability of the code enables incorporation of plasma actuators into the computations with relative ease 
since electrodes can be defined as individual solid blocks. The domain can be divided into two separate 
computational domains: one for the air side and the other for the dielectric wall. The GHOST code has been 
previously validated against a wide range of test cases and flow conditions and has been used extensively in several 
low pressure turbine related publications. 30 ' 33 


III. Results and Discussion 

A. Model Validation and Calibration using Quiescent Flow Experiments 

In order to calibrate the parameters 
Debye length, X d and the maximum 

charge density on the wall, , appearing 
in the model we have chosen the quiescent 
flow experiments (Ref 34) conducted using 
a single pair of electrodes to better isolate 
the effects of the actuator on the 
surrounding air. The details of the actuator 
geometry and experimental set up are 
shown in Figure 6. The actuator consists 
of two 1cm wide conductive copper strips 
as electrodes which are separated by 
Kapton dielectric with e r value of 2. 7. 

In the experiments the lower electrode 
is grounded and plasma region is generated 
using a square wave with frequency of, d) = 4.5 kHz and amplitude of <p maK = 5 kV. The experimental data include 
velocity field measurements obtained using multiple PIV systems. Further details of the experiments and the 
measurement techniques used are given in Ref. 34. It should be noted that the experimental data are not intended to 
be used to calibrate the model and the experimental errors have not been quantified. Therefore, we shall only use 
these data for purpose of demonstration and a more comprehensive refinement of the model will follow in future 
work. 


t e - 0.102 mm x e = 0.5 mm 
L e =10mm y e = 0.127 mm 



Figure 6. Experimental set up and details of the plasma actuator. 
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In the computations our aim is to match the maximum velocity observed in the experiments as well as the 
experimentally observed flow pattern shown in Figure 7a. From the experiment it is observed that the flow is drawn 
into the surface region above the embedded electrode by the plasma induced body force. This results in a jet issuing 
to the right of the actuator with a maximum velocity of approximately 1 m/s. 



(a) Experiment 34 



Figure 7. Comparison of experimental 34 and computed streamlines for the plasma actuator in quiescent flow. 
The actuator interface is located at the 2cm tick. 


Based on the flow pattern and maximum velocity criteria, a choice of the parameter values of A, =0.00 lm and 


p™ x =0.0008C/m 3 seems to agree well with the experiment. The streamlines obtained from the computation using 


The computed electric potential distribution in the 
vicinity of the electrodes obtained from equation (7) is 
shown in Figure 8 along with the streamlines of the 
actuator induced flow. The computed electric 
potentials show that the strongest electric potential 
variation, or the electric field, is in the region between 
the two electrodes. This is also the region where the 
strongest concentration of plasma is observed in the 
experiments. The streamlines indicate that the flow is 
pulled from above into this region and jetted to the 
right direction. 


Figure 8. Computed electric potential contours and 
streamlines in the vicinity of the electrodes. 
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these values are shown in Figure (7b). Although not an exact match, the computed flow pattern compares favorably 
with the experimental flowfield shown in Figure (7a). The computed maximum velocity is lm/s also matching the 
experimental value. It should be noted that the boundary layer obtained by the computation appears to be thinner 
than the experimental data. Although the thickness 
can be adjusted by increasing the value of the Debye 
length, X d , we did not attempt to do so because the 


experimental data may be contaminated by 3-D effects 
caused by the experimental setup. A refinement of the 
model will be performed when more accurate data are 
available. 
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Figure 9. Computed charge density contours in the 
vicinity of the electrodes. 


The charge density distribution obtained from 
equation (19) is shown in Figure 9. The charge 
density has its maximum value- in the region over the 
comer of the embedded electrode as specified by the 
half Gaussian distribution. A thin layer of plasma is 
formed on the surface above the embedded electrode 
and gradually creeps into the gap staggered between 
the two electrodes as shown in Figure 9. This 
predicted charge density distribution resembles the 
plasma profile observed in the experiments. 

The contours of magnitude of the computed 
body force vectors are shown in Figure 10. The 
highest body force is concentrated in the region over 
the upper left comer of the embedded electrode, 
where the concentration of charge density and the 
electric field strength are at maximum, consistent 
with the location where the strongest plasma 
generation is observed experimentally. 


These results indicate that the model can mimic 
the experimentally observed effects caused by 
plasma actuators. It must be noted that, in this study 
we have limited the calibration of the model to a 
single experiment in order to provide a proof of 
concept for the model and the modeling approach. 
The results illustrate that the model and the 
approach are promising in the computation of 
plasma flow control applications. Flowever, the 
model has to be further tested and the model 
parameters need to be further calibrated using more 
reliable experimental data. 


B. Flow Separation Control in a Low Pressure 
Turbine Cascade 

In order to illustrate the models applicability to 
plasma flow control applications a test case 
involving a low pressure turbine cascade is 
computed. We have selected an intermediate Reynolds number case from the Pratt and Whitney Pak-B blade 
experiments of Huang et al. 9 * . We chose this case for demonstration purposes because it has a large flow separation 
region and the experimental data taken includes velocity profiles at various stations on the suction side of the blade 
and the pressure coefficient distribution. Further comparisons are available in Jacob et ai. 16 The flow conditions are 
Re= 50,000 based on cascade inlet velocity and axial chord length and the freestream turbulence intensity is 0.08%. 
The experimental data for the base flow without flow control was used in this study. The flow separation zone is 
characterized by the plateau in the Cp distribution on the suction side of the blade. It is observed experimentally that 
the separation for the no-actuator flow extends from x/C s ~0.7 to 0.95, as shown in Figure 11. Figure 11 also 
shows the comparison of the C p distribution on the surface of the blade from the computation. The results show the 
separation starts at x / C x =0.7 but ends slightly earlier than the experimental location. 



Figure 10. Computed body force magnitude contours in 
the vicinity of the electrodes. 
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Figure 11. Comparison of pressure coefficient 
distribution for Pak-B Blade, Re=50,000, 
FSTI=0.08%. 


The same actuator configuration as discussed in 
section IIIA is introduced on the suction surface of the 
blade with actuator interface located at 65% axial 
chord. Note that the actuator that was used in the 
simulation is not identical to the actuator used in 
Huang et al. 9 A sine wave voltage with a frequency of 
5kHz, and an amplitude of 5kV is applied in the 
current demonstration. As can be seen from Figure 11, 
the use of the plasma actuator has caused the 
separation to reduce in size. The comparison of 
velocity profiles at various stations on the suction 
surface of the blade is shown in Figure 12. The figure 
shows that the size of the separation is reduced to end 
at ,r/C r =0.87. Comparisons with experimental data 
of the same geometry using a plasma actuator show 
similar results, in that separation is not completely 
eliminated and an inflection point is still present, but 
losses are significantly reduced downstream of the 
actuator location. 16 The ability of the model to 
simulate a control of flow separation for a low 
pressure turbine blade is thus demonstrated. 


LV. Concluding Remarks 

The current paper presents a new modeling approach to simulate flow control using plasma actuators. The 
model solves two additional equations representing electrical field and charge density of the ionized air. These two 
equations need to be solved only once before the computation. A body force vector can be deduced from the 
solutions of these two quantities and it can then be incorporated into Navier-Stokes equations to account for the 
plasma-actuator effects. The model is calibrated against a simple plasma-actuator-driven flow in a quiescent 
environment. The model is then used to simulate an actual low-pressure turbine flow to achieve reduction of flow 
separation. 
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